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Abstract 

A detailed understanding of all effects and influences on higher-order correlations is essential. At 
low charged multiplicity, the effect of a nonpoissonian multiplicity distribution can significantly dis- 
tort correlations. Evidently, the reference samples with respect to which correlations are measured 
should yield a null result in the absence of correlations. We show how the careful specification of 
desired properties necessarily leads to an average-of-multinomials reference sample. The resulting 
internal cumulants and their averaging over several multiplicities fulfil all requirements of correctly 
taking into account nonpoissonian multiplicity distributions as well as yielding a null result for 
uncorrelated fixed- iV samples. Various correction factors are shown to be approximations at best. 
Careful rederivation of statistical variances and covariances within the frequentist approach yields 
errors for cumulants that differ from those used so far. We finally briefly discuss the implementation 
of the analysis through a multiple event buffer algorithm. 
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1 Introduction and motivation 



The understanding of hadronic colUsions is now considered an essential baseline for ultrarelativistic 
heavy-ion collisions. Given the correspondingly low final-state multiplicities, there are significant de- 
viations, even for inclusive samples, from assumptions commonly made both in the general theory 
and in the definition of experimentally measured quantities such as a nongaussian shape of the corre- 
lation function and nonpoissonian multiplicity distributions. Constraints such as energy-momentum 
conservation [H [2] would also play a role in at least some regions of phase space. Multiplicity-class 
and fixed-multiplicity analysis differ increasingly from poissonian and inclusive distributions, and with 
the good statistics now available, measurements have become accurate enough to require proper un- 
derstanding and treatment of these assumptions and deviations, which play an ever larger role with 
increasing order of correlation. 



1.1 Correlations as a function of charged multiplicity 

There are a number of reasons to study correlations at fixed charged multiplicity or, if necessary, 
charged-multiplicity classes. Firstly, the physics of multiparticle correlations will evidently change with 
N, and indeed the multiplicity dependence of various quantities such as the intercept parameter and 
radii associated with Gaussian parametrisations is under constant scrutiny [3l [31 [5l [6l [71 [Hj [9l [TOl [TT]. 
Measurement of many observables as a function of multiplicity class, regarded a proxy for centrality 
dependence, has been routine for years. Corresponding theoretical considerations, e.g. in the quantum 
optical approach go back a long time [12]. Secondly, correlations for fixed- are the building blocks 
which are combined into multiplicity-class- and inclusive correlations [13j. 

However, such fixed- correlations have been beset by an inconsistency in that they are nonzero 
even when the underlying sample is uncorrelated and do not integrate to zero either. This has been 
recognised from the start [H], and various attempts have been made to fix the problem. 

Combining events from several fixed- A/^ subsamples into multiplicity classes does not solve these 
problems. To quote an early reference [15]: Averaging over multiplicities inextricably mixes the prop- 
erties of the correlation mechanism with those of the multiplicity distribution. Instead, the study of 
correlations at fixed multiplicities allows one to separate both effects and to investigate the behaviour of 
correlation functions as a function of multiplicity. Under the somewhat inappropriate name of "Long- 
Range Short-Range correlations" |141 [T6] , an attempt was made to separate these multiplicity-mixing 
correlations from the fixed- A/^ correlations, but the inconsistencies inherent in the underlying fixed- A^ 
correlations were not addressed. Building on Ref. [17], we propose doing so now. 



1.2 Cumulants in multiparticle physics 

Multiparticle cumulants have entered the mainstream of analysis, as shown by the following incomplete 
list of topics. In principle, the considerations presented in this paper would apply to any and all such 
cumulants to the degree that their reference distribution deviates from a poisson process or that the 
type of particle kept fixed differs from the particle being analysed. 

Integrated cumulants of multiplicity distributions have a long history in multiparticle physics 
|18j . Second-order differential cumulants, normally termed "correlation functions", have likewise been 
ubiquitous for decades [6] both in charged-particles correlations [l^ and in femtoscopy since they 
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provide information on spacetime characteristics of the emitting sources, most recently at the LHC 
[9l [TOl [T9] . Differential three-particle cumulants generically measure asymmetries in source geometry 
and exchange amplitude phases ^20j . They also provide consistency checks [2]J and a tool to disentangle 
the coherence parameter from other effects \22\ I23j . Three-particle cumulants are also sensitive to 
differences between longitudinal and transverse correlation lengths in the Lund model |24) . Inclusive 
three-particle cumulants have been measured, albeit with different methodologies, in, for example, 
hadronic [2"51 Eg [271 [2S] , leptonic [Ml EHl ED 132] and nuclear coUisions [Ml [Ml ESI ES] • They play a 
central role in direct QCD-based calculations [371 EEl [39] and in some recent theory and experiment 
of azimuthal and jet-like correlations [4U\ WI\ WI[ \^[ I44| . Net-charge and other charge combinations 
are considered probes of the QCD phase diagram [451 1^61 147] . Cumulants of order 4 or higher are, 
of course, increasingly difficult to measure and so early investigations were largely confined to their 
scale dependence UHl 091 [501 [51] . The large event samples now available have, however, made feasible 
measurements of fourth- and higher-order cumulants in other variables as proposed in |52l \53\ [Ml 
[55l 156] as, for example, recently measured by ALICE [57]. Reviews of femtoscopy theory range from 
[551 l59l l60j to more recent ones such as [7j . 

I. 3 Outline of this paper 

It has long been obvious that the root cause of the problems and inconsistencies set out in Section 

II. 11 was the reference sample [61j. Insofar as cumulants are concerned, the solution was outlined in 
Ref. [T7| as a subtraction of the reference sample cumulant from the measured one; important pieces 
of the puzzle were, however, still missing at that stage. In this paper, we clarify and extend the basic 
concept of internal cumulants and consider in detail the case of second- and third-order differential 
cumulants in the invariant Q = \/—{pi — ^2)^ for fixed charged multiplicity A^. The method may be 
implemented for other variables without much fuss. 

A second cornerstone of the present paper is the recognition that the n particles which enter a 
correlation analysis are usually only a subset of the A^ charged pions. While in the case of charged- 
particle correlations all A^ particles are used in the analysis, Bose-Einstein correlations, for example, 
would use only the n = positive pions (and, in a separate analysis, only the n_ = N—n^ negatives). 
In addition, there may be reasons to restrict the analysis itself to subregions of the total acceptance 
Q in which A^ was measured, as exemplified in this paper by restriction to a "good azimuthal region" 
subinterval around the beam axis, A C [0,27r], in which detection efficiency is high. A can, however, 
be reinterpreted generically as any restriction in momentum space compared to and/or as a selection 
such as charge or particle species. Even when setting A = i.e. doing the femtoscopy analysis in 
the full acceptance n still does not equal A^ but fluctuates around A^/2. The trivial observation that 
n ^ N fundamentally changes the analysis: identical-particle correlations at fixed N and charged- 
particle correlations at fixed N require different definitions. 

As we shall show, ad hoc prescriptions such as simply inserting prefactors or implementing event 
mixing using only events of the same A^ do alleviate the effect of the overall nonpoissonian multiplicity 
distribution in part but fail to remove them completely. The same issues will, of course, arise in 
any other correlation type of, for example, nonidentical particles or net charge correlations. The 
formalism set out here can be easily extended to such cases. A refined version of the abovementioned 
Long-Range-Short-Range method, which we term "Averaged-Internal" cumulants, will be presented 
in Section [5] Along the way, we document in Section [2] extended versions of the particle counters 
|62[ [55] which we need as the basis for correlation studies and in Section [3] demonstrate from first 
principles that statistical errors for cumulants used so far have captured only some of the terms and 
with partly incorrect prefactors. Section [U] outlines the implementation of event mixing for fixed- A^ 
analysis. While experimental results will be published elsewhere, preliminary results in Figs. 2 and 3 
show that, in third and even in second order, corrections due to proper treatment of fixed- A^ reference 
samples can be large. 
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2 Raw data, counters and densities 



2.1 Raw data 

The starting point for experimental correlation analysis is the inclusive sample S, made up of £ 
events a = Each event consists of a varying number of final-state elementary particles 

and photons; for our purposes, we consider only the N{a) charged pions of event a in fi, the 
maximal acceptance region used. Each pion i = l,...,N(a) is characterised by a data vector 
{P^,sf,e'^) containing its measured information, including the three components of its momentum, 
-P" = iPixiPiyPtz) (yf ^ 4'i iPti)^ while its discrete attributes such as mass, charge etc. are captured 
in a data vector sf of discrete values; for the moment, we shall consider only the charge, sf — )• cf. 
From the sample's raw data, we can immediately find derived quantities such as the total charged 
multiplicity N{a), the total transverse energy etc., and such derived quantities are hence considered 
part of the raw data. The list of particle attributes should be augmented by an error vector ef con- 
taining the measurement errors for each track, but we shall not consider detector resolution errors 
here. In summary, the inclusive data sample is fully described in terms of 5 = {P", s"'}^^^ consisting 
of lists of vectors in continuous and discrete spaces 

P*^ = {Pf , . . . , P^(^^} = {s", • • • 5 •SAr(a)}- (1) 

2.2 Data after conditioning and cuts 

For a particular analysis, the inclusive sample is invariably subdivided and modified through "con- 
ditioning", the statistics terminology for semi- inclusive or triggered analysis: From the total sample 
of events, a subsample is selected according to some restriction or precondition. In our case, this 
conditioning proceeds in the following steps: 

• Conditioning into fixed-N subsamples: For the fixed-multiplicity analyses that form the subject 
of this paper, S is subdivided into a set of fixed- A'' subsamples S^, each of which contains only 
events a whose measured multiplicity N{a) is equal to the constant N characterising 5jv = 
{P", I 5{N, N(a))}, N = 0,1,2, .. . We use the vertical bar | here and everywhere in the usual 
sense of "conditioning" whereby the events in sample 5jv must satisfy the condition that their 
charged multiplicity must equal the specified constant N, denoted in this case by the Kronecker 
delta d{N, N{a)). Quantities to the right of the vertical bar are generally considered known and 
fixed, while quantities left of the bar are variable or unknown. The number of events in 5jv 
equals the (5-restricted sum over the inclusive sample, 

S 00 
£^ = Y,m,N{a)), ^£^ = £. (2) 

a=l Af=0 

The usual multiplicity distribution is the list of relative frequencies 

N=0 

While desirable, it is not easy to measure the total multiplicity of final-state charged pions, 
a quantity which approximately tracks the variation in the physics. Choosing charged pions 
measured within the maximal detector acceptance N = N{Q) as marker is in any case only an 
approximation because it excludes charged particles outside the primary cuts and also ignores 

^ These event ratios are not probabilities in the strict sense; in the frequentist definition of probability, the two are 
equal only in the limit £ — >■ 00. We therefore avoid the use of the symbol Pn for such and similar data ratios. 
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final-state particles other than charged pions. Nevertheless, we expect to be a reasonable 
measure of the multiplicity dependence of the physics. Alternatively, the multiplicity density in 
pseudorapidity at central rapidities dN/drj can be used as a model-dependent proxy for N. 

• Azimuthal cut: While N(a) is the charged multiplicity measured in 17, there is no a priori 
reason why the correlation analysis itself may not be conducted within a restricted part A G 
of momentum space within which the actual analysis is done. In the case of the UAl detector 
from which the data used in the examples below was drawn, A refers to azimuthal regions within 
which measurement efficiency was high, and pions found in the low-efficiency azimuthal regions 
were excluded. Correspondingly, the multiplicity n{a) which enters the correlation analysis itself 
differs from N{a) and will, for a given fixed fluctuate with relative frequency 

T^nN = -7; — for each fixed A'" = 1, 2, ... , (4) 

where £nN is the number of events for which N{a) = N and n(a) = n. The outcome space 
for n(a) will depend on its definition; in the present case where only positive (or only negative) 
pions within A are used in the analysis, it will be [0, 1, . . . , A^] so that the relative frequency is 
normalised by 

N 

J^7^„^ = l VAT =1,2,.... (5) 

n=0 

With approximate charge conservation n-|- ~ n_, we expect the fixed- A^ average for positive (or 
negative) pions in A to hover around 

^ ^ A^ volume of A 
2 volume of 17 

An example of the resulting relative frequencies (conditional normalised multiplicity distribu- 
tions) is shown in Fig. 1. Since n < N, these conditional multiplicity distributions are almost 
always subpoissonian, i.e. narrower than a Poisson distribution with the same {n)^ would be. 

• Generalisation: While in this paper the analysis will be carried out for the n(a) positive pions of 
event a falling into A, the same formalism obviously applies to negative pions and may equally 
refer to any other particles such as kaons, baryons, photons etc in any combination in dependence 
on A^. There is no a priori connection between the definitions of A^ and n. 

• Identical-particle vs multi-species analysis: While we do not develop the formalism for correla- 
tions between two or three particles of different species or charge, the methodology developed 
here can be easily modified to deal with such cases. For example, positive-negative pion com- 
binations and "charge balance correlations" [65] can be handled by inserting delta functions 
(5(c, c°), where c is the desired charge and cf the measured charge of track i in event a, into the 
definitions of the counters in Section 12.31 

2.3 Counters and densities for fixed A^ 

This section is based on an old formalism |62 1I631[52] which must, however, be updated to accommodate 
the issues being considered here. The basic building block of correlation analysis is the counter; it is 
a particular projection of the raw data particular suited to the construction of histograms. Eventwise 
counters p for a given event a are averaged to give sample counters p. 



5 




Figure 1: Conditional normalised multiplicity distribution (relative frequency) of the number of positive pions 
n in restricted azimuthal region A = {20° < < 160°} for respectively fixed charged multiplicity iV = 3, 5, 
10, 20, 30, 40, 50, for the UAl dataset used in Ref. [M]- In accordance with Eq. ©, (n)^ ~ 0.39iV. 



We take the simple case where event a contains N(a) tracks with three-momenta = {Pf", ■ ■ ■ , P^^^-^}, 
no discrete attributes s and no further cuts or selection. For each point pi in momentum space, only 
that particle i (if any) whose momentum Pf happens to coincide with p is to be counted, 

N{a) 

p{p^\p'^) = Y,5{p,-pn. (7) 

1=1 

Such counters always appear under an integral over some region of the -P-space, so that the delta 
functions fulfil the purpose of counting those particles falling within that region. Alternatively, one 
can consider the delta functions here and below to represent small nonoverlapping intervals around 
the specified momenta. The integral over the full momentum space yields 

/ dpi p{pi I P") = N{a) (8) 
Jq. 

while an integral over some subspace or bin C $1 will yield the number of particles of event a in 
bin Qh. The second-order eventwise counter for event a iqfl 

Nia) 
i^j=l 

with the inequality i ^ j ensuring that a single particle is not counted as a "pair". The counter 
integrates to 

[ dpidp2p{pi,P2\Pn = N{af (10) 
Jn 



^Pairs are ordered, i.e. a particular pair is counted twice. Unordered pair counting is possible but unnecessarily 
complicates sum limits. 



6 



using the falling factorial notation 

N^ = N{N -1)---{N -r + 1) r = l,2,... (11) 

as contrasted to the rising factorial (Pochhammer symbol) A^*" = N{N + 1) • • • {N + r — 1). The 
single-particle counter is a projection of p{pi,P2 \ -P") because 

/ dp2 p{Pi,P2 I PI = [N{a) - 1] p{pi I P-^). (12) 
Jn 

The most general eventwise counter which enters the exclusive cross section for events with charged 
multiplicity N{a) 

N{a) N{a) 

p{Pl,P2,...,PNia)\Pl= ^ U^^Pd-K) (13) 

jl7^i27^---7^«jV(a) = l d=l 

fully describes the event, including any and all correlations between its particles. It integrates to the 
factorial of the event multiplicity 

/ dpidp2 ■■■ dpiy(^a) p{Pl^P2,---,PNia)\P'') = N{a)\ (14) 

Jn 

and contains all counters of lower order by projection. An rth-order counter p{pi,P2, ■ ■ ,iPr \ P"') is 
zero whenever there are more observation points than particles being observed, r > A^(o)o 

To distinguish eventwise counters for nonfixed N from eventwise counters for fixed N, we define 
the separate eventwise counter for fixed N by specifying an additional Kronecker delta, 

N 

p{pi,...,Pr\P'',N) = 6{N,N{a)) Yl Il^(Pd-P^") r = l,2,...,N. (15) 

While the counters and densities defined above and below are clearly frame-dependent, it is easy to 
define corresponding Lorentz-invariant versions by supplementing each delta function in 3-momentum 
with the corresponding energy; thus Eq. ([7]) would become, for example, 

N(a) 

p{p^\P-)=Y,E{pi)5{p^-Pt) (16) 



with E{jpi) = El = \/ Pi + rrfl the on-shell energy, and in general 

N{a) N{a) 

p{puP2,...,Pr\Pl= Yl llEd6{pd-P0 r = l,...,N{a), (17) 

N 

p{pi,...,Pr\P'',N) = 6{N,N{a)) Y llEdd{pd-Pt,) r = l,...,iV, (18) 

which are manifestly invariant. Because such counters and densities are, however, always integrated 
over some by Y\^{dpd/ Ed), the additional factors E^ always cancel and play no role on this level 
of analysis and will be ignored for the time being. The bin boundaries of il;, do, however, remain 
frame-dependent. 



counter of order N{a) can be made to behave like one of order N{a) + 1 by defining an additional dummy data 
point PN{a)+i which lies outside the normal domain f2, but we shall not pursue this here.. 
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Charge-, spin- or species-specific counters are defined in the same way, i.e. by supplying appropriate 
Kronecker deltas to the counters; for example the particle counter for pions with charge ci at pi for 
fixed N is 

Nia) 

pipi\ci,P'',N) = 6iN,Nia)) ^ 5(ci, <5(pi - i^") (19) 

i=l 

while the two-particle counter for charge combination (01,02) at momenta (pi,P2) for any is, for 
example, 

Nia) 

P{PI,P2 I Cl, C2, P-^) = '^(Cl' -^(Pl - Pi) ^^(^2, C,^) 5{P2 - P/). (20) 



In contrast to Eq. ()19p . charge counters rather than particle counters would be 

N{a) 



,{pi \ci,P-)=Y, ci <5(ci, S{pi - Pt) (21) 



1=1 



so that Pc{pi I +1,P") + Pc{Pi I — IjP") represents the net charge of event a at pi. The two-particle 
counter for charges {01,02) at momenta (pi,P2) is 

Nia) 

p'^'^' = p,{pi,p2\ou02,Pn = Y ci5ioi,o'})5{pi- Pno25{o2,o';)5{p2- P,^), (22) 

and "charge flow" correlations can be constructed from this (for rapidities {y, y') in the case of Ref. 
[66] ) such as 

Hy, y') = -(E,,^, ct 0- 5{y - Y-) 5{y' - Y-)) (23) 

which can be expressed as <I> = (p^ + p^ — — p ) and the related "charge balance functions" 
described in e.g. [45j. 

Returning to the fixed- case: eventwise counters will usually be combined with similar events to 
form event averages. The simplest average is the fixed- density for the subsample of fixed 5jv, 

p(^>l|5^.) = -3-5;p(Pl|P^7V), (24) 

a 

with the Kronecker delta in ()15p ensuring that only events in 5jv are considered, so we need not further 
specify the individual terms or limits of the a-sum. Using ([2]), it is immediately clear that 

/ dpi p{pi \S^) = N (25) 

compared to the integral over the corresponding eventwise counter 

[ dpp{p\P'',N) = NS{N,N{a)) (26) 
Jn 

and to the integral ([5]); similarly 

/ dpi--- dpr p{pi, ...,Pr\ S^) = N{N-l) - - - {N-r + 1) = A^^. (27) 
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The inclusive averaged density p{pi \ S) is the weighted average over all N of the fixed- averages, 

oo 

p(pilcS) = ^7^^p(pl|5^). (28) 

N=l 

Using ([3]l. (fT5]) and ([2l|) . this can be written as 

a i 

and for general r = 1,2,... 

oo \ ^ 

p{pi,...,Pr\S)=Y,T^Np{Pl,---,Pr\S^) = -Y, Yl Il^(Pd-K)^ (30) 

N=r a d=l 

keeping in mind that p{pi, . . . ,Pr\ 5jv) will be zero whenever N{a) < r. The integral of any rth order 
inclusive averaged density is the rth-order factorial moment of the multiplicity distribution, 

„ oo 

/ dpi--- dpr p{pi, ...,Pr\S)=J2'^NN'-= (iV-) (31) 

with simple angle brackets denoting inclusive averaging. 

The averaged counters are of course directly related to the traditional definitions in terms of 
cross sections. If C is the integrated luminosity of incoming particles, the topological cross section is 
cTjv = Sn/J--, the inelastic cross section is ctj = £/C and the inclusive cross section is crinci = J2n = 
(A^) a I while the relative frequency (multiplicity distribution) can be written as usual as TZ^ = ctm/cti- 
The relation between the differential cross sections and our counters is 

Pmcl(Pl,---,Piv) = p{pi,. . . ,Pn\S) = —-^ ^mcl_ ^^2) 

0-j dpi--- dpN 

, I > 1 d ^^Cxcl fe\f\\ 

PPi,...,Pjv = — — , (33) 

cjjv dpi--- dpN 

and so as usual inclusive and exclusive densities are related by [13] 

p{pi, ...,Pn\S) =Y^ / p{pi, . . . ,Pjv I 5jv) dpr+i - - - dpj^, (34) 

N=r ^ 

while the semi-inclusive cross sections and counters follow by the usual projections. 



2.4 Counters and densities for fixed {N,n) 

Our choice of a basic counter is motivated by the experimental situation set out in Section [1} we wish 
to work in event subsamples of fixed total charged multiplicity N{a) in the entire momentum space 0, 
but do the differential correlation analysis using only those pions n which fall into the restricted space 
A and of a particular charge +1 or —1. This requires the use of "subsubsamples" for which both N 
and n are kept fixed, 

•SnN = {P°', with a constrained by 5(n,?i+(a)) 5{N, N{a))~\ (35) 
with n+(a) the number of positive pions of event a in A, and eventwise subsubsample counters 

n 

p{pi,...,Pr\n,N,P-) = 6{n,n^{a))6{N,N{a)) J] S{pi - Pt,) - - - S{pr - Pt,.) . (36) 

^lJ^■■■J^^r='^ 
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As in Eq. the number of events in a subsubsample £nN = '^(^i ^+(^)) '^(-^' ^('^)) writers the 
relevant event averages 

p(pi,...,Pr\SnN) = -^^p{pi,...,Pr\n,N,P'') (37) 

a 

where once again the double Kronecker deltas in p6p ensure selection of events in SnN only. Integrals 
of the counters over the good-azimuth region A yield, for the eventwise and 5„jv-averaged counters, 

/ dpi---dpr p{pi,...,Pr\n,N,P'') = n^6{n,n+{a)) 6{N,N{a)) (38) 

Ja 

dpi--- dpr p{pi , . . . , Pr I 5„jv) = n-. (39) 

Bearing in mind that observation points pi,p2, ■ ■ ■ refer to positive pions in A only, the event-averaged 
counters for fixed N but any n are given by the average weighted in terms of the relative frequency 
"^71 JV = Sun/Sni 

N 

p{pi, ...,Pr\ Sn) = ^ TZnN P(P1, • • • ,Pr I 5„iv) = (p(Pl, • • • , Pr | ^nAr))^ (40) 
n=r 

for r = 1, 2, 3, . . ., which integrate to 

dpi--- dpr p{pi, ...,Pr\S,^) = (n-)^. (41) 

3 Construction of correlation quantities 
3.1 Criteria 

Correlation measurements of any sort are only meaningful if a reference baseline signifying "indepen- 
dence" or "lack of correlation" is defined quantitatively; indeed, many different kinds of correlations 
may be defined and measured on the same data, depending on which particular physical and mathe- 
matical scenario is considered to be known or trivial and taken to be the baseline [HT]. In our case, 
we require the reference distribution to have the following properties: 

1. The number of charged pions in all phase space N is an important parameter as a measure of 
possibly different physics, but only the n positive pions in A are to be considered in the differential 
analysis. 

2. For a given {N, n), the momenta of the reference density p'^'^^{pi, . . . ,pr\ Sun) should be mutually 
independent for any order 1 < r < n. This and the previous requirement imply that the reference 
should be a n-multinomial distributed over continuous momentum space; see Section 13.2.11 

3. Given fixed N, the reference density p^'^^{pi, . . . ,Pr \ S^) must reproduce the n -multiplicity struc- 
ture of the subsubsamples Sun o-s embodied in IZun- As set out further in Section 13.2.21 this 
translates into an average of multinomials, 

N 

p'^'ipi, ...,Pr\ Sj,) = Y,'^nN P^'^'^PU - - - ,Pr | «,5„iv) = (p'""'*(Pl, • • • , Pr | CX.Sn^))^. (42) 

n=r 

4. The reference density should reproduce the measured one-particle density in momentum space. 
This can in principle be satisfied by three different expressions for the multinomial's parameters 
a: see Section [3.2.31 
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5. Measures of correlation must reduce to zero even on a differential basis whenever the data is, in 
fact, uncorrelated. While this may seem self-evident, this requirement is often ignored or not 
satisfied in the literature. We address the resulting proper baseline through the use of internal 
cumulants in Section [3131 

6. The measure of correlation should he insensitive to the one-particle distribution. This is addressed 
as usual by normalisation; see Section I3.3[ 



3.2 The reference distribution 

3.2.1 Multinomials in discrete and continuous spaces 



Before Eq. (p2|) can be developed further, it is necessary to take a detour into discrete outcome spaces 
before tackling the continuous outcome space defined by p and P"^. The reason is that multinomial 
distributions for continuous arguments p can be written only as a limit of the discrete precursor. 

Let there be bins fib, 6 = with the corresponding set of Bernoulli probabilities a = 

{a(6)}^-i^ of a single particle falling into bin fi^, normalised by '^i,ci{b) = 1. Independent tossing of 
n particles into these bins results in the multinomial for the bin counts n = {j^fe}^^. 



B 



p{n I a, n) = n\ 



a{b) 



rib 



6=1 



with normalisation 



p{n\cx,n) = 1, 

r/(n) 

where the sum must be taken over the "universal set" 

U{n) = {n 1 > 0; Ef, = n}. 



(43) 



(44) 



(45) 



The multivariate factorial moment generating function (FMGF) for this multinomial for the set of 
source parameters A = {A(6)}^^ can be solved in closed form. 



Q^-'\\\a,n) = ^p{n\a,n)ll{l-Xib)r = l-^ma{b) 



(46) 



The FMGF Q{X) can generally be used to find multivariate factorial moments p{bi^ ,bi^, - ■ ■ ,bi^) and 
factorial cumulants K^bi^^ , bi^ , ■ ■ ■ ,bi^.) for any selection of bins {bi^^ , bi^ , ■ ■ ■ ,bi^.) G (1,... ,B), including 
repeated indices, by differentiation 



2 ' ' "l-r 



l^ibii 1 bi 



dx{bi,) dx{bi,) ■■■dx{bi^) 

{-ly aMnQ(A) 



12 1 1 "Ir 



dx{bi,)dx{h 



«2 ; 



5A(6iJ 



A=0 



A=0 



For the multinomial case (|46|) . the factorial moments and cumulants are therefore 

r 

p'^'''\b,„b,„--- ,bijcx,n)=n^ lla{id) ^ r < n, 

d=l 

r 

Ac-"i*(6,, , • • • , 6iJ a, n) = (-1)^"^ (r - 1)! • n "(^d)- 



(47) 
(48) 

(49) 
(50) 



d=l 
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The multinomial for variable p in continuous outcome space M is derived by keeping n constant 
while taking the limit B ^ oo with bin sizes tending to zero and changing to a Bernoulli probability 
density a{b) — >• dpa{p) normalised by J^dpa{p) = 1. The result is the point process where the 
probability for the count n(p) in the infinitesimal "bin" around any p to be larger than 1 becomes 
negligible, i.e. we have at most one particle at a given p. While the multinomial probability itself can 
be written only as a limit, the FMGF can be written analytically as the functional [67] 



Q^''''[X{p)\aip),n] 



1 



A 



dp\{p)a{p) 



Factorial moments and factorial cumulants are found generically from functional derivatives [13 

(-1)^ S'Q[\{p)] 



l^iPil ; Pi2 ) ■ ■ ■ ; Pir] 



{-If <5nnQ[A(p)] 



A(p)=0 
A(p)=o' 



(51) 

(52) 
(53) 



which for the multinomial Q[\{p)] = Q™^^^'[A(p) | a{p),n\ of ([5T]) yield 

r 

/0™"^*(Pii,Pi2, • • • ^Pir I oc{p),n) =n-W a{pij 1 < r < n, 

k=l 

r 

«""'*(Pn,Pi2r ■■ ,P^r\ = (-1)^-' (r - 1)! • n n 



(54) 
(55) 



k=l 



3.2.2 Multinomial reference for fixed N 



Applying the above general case to our reference distribution (|42|) . we must rewrite Eq. (|51|) to make 
provision for the fact that a may in general depend not only on A'^ but also on n, 



Q-"'*[A(p)|a(p|5„^)] 



1 - y dpX{p)a{p\S„ 



(56) 



Inserting ()56p into ()42p . we find the FMGF for the reference distribution of subsample to be 



Q'"'[A(p) I «(P), 5^)] = Q'^'^'^Kp) I a{p I Sn^)] 

n 

Using ()54p . the reference factorial moments are therefore, 



I - j dpX{p)a{p\Sn 



■ (57) 



P^''^{Pl,P2,--- ,Pt\Si, 



n- 



Y\ a{Pk I Sun) 



\/ r < n 



(58) 



\ fc=l In 

with corresponding expressions for the reference factorial cumulants. 



3.2.3 Reproducing the one-particle distribution 

The set of functions a{p \ Sun) are as yet undetermined, apart from the general constraints a{p \ Sun) > 
and dpa{p \ SnN) = 1. In multinomials of all kinds, the Bernoulli probabilities cx are fixed param- 
eters and therefore are the conveyers of whatever remains constant in the outcomes while the detailed 
outcomes fluctuate as statistical outcomes do. The "field" a(p|5„iv) can and must therefore be seen 
as the quantity encompassing the "physics" of the one-particle distributions, which, in the absence of 
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additional external information, is embodied by our experimental data sample: the experimental den- 
sities p{pi, ■ ■ ■ ,Pn I <Sn) "are" the physics, including all correlations, and their first-order projections 
p{pi |5jv) "are" the one-particle physics. The question immediately arises whether a{p\SnN) should 
be fixed by p{p \ Snn) or the n-average p{p \ S^) = {pip \ SnN))i^- Three possible choices come to mind: 



1. It is tempting to define it in terms of the density for each subsubsample 5, 



nNi 



a{p\S^,) = M^ V(iV,n), (59) 

which is correctly normalised since J dp p{p \ Sun) = n. As this choice would attribute physical 
significance to n, it would be appropriate whenever n is associated with additionally measured 
experimental information. If, however, n fluctuates randomly from event to event based in part 
on unmeasured or unmeasurable properties such as an event's azimuthal orientation, use of ()59p 
makes no sense. 

If n is deemed physically relevant, correlations in terms of p{pi, . . . ,Pr\ 5„jv) of Eq. (f37|l may be 
feasible, conditional on the availability of a sufficient number of events Sun- Where sample sizes 
do not permit this, one could nevertheless attempt to measure what have historically been termed 
"short-range correlations" but in this case not in the traditional sense of fixed-N correlations 
versus inclusive ones, but rather of fixed-n- fixed- correlations versus fluctuating-n- fixed- A^ 
correlations. See Section [3^61 



2. A second choice 



a(p|5..) = /^(^^\ (60) 

\ ^ /n 

would be properly normalised but fails to satisfy the crucial relations (|78p ~ (j80p below and is 
hence discarded. 

While remaining open-minded towards Choice 1, we therefore choose the third possibility, the 
ratio of the average density divided by the average, all for fixed A^, 

which would be appropriate for samples where £nN is too small or physical significance can be 
attributed only to A^ but not to n. According to (j4ip . it is also correctly normalised and ensures 
that the Bernoulli parameters are the same for all events in S^, independent of n. Substituting 
this into Eq. (j58p . the differential reference factorial moments orders become 

p"'(pi,P2, • • • ,Pr 1 5^) = 7^ n p^Pd I = n p^P'i I '^^) (62) 

^"^^^ d=i d=i 

where we identify the prefactor as the normalised factorial moments of the n-multiplicity distri- 
bution for given A^, 

(n-) 

FrN = , , r 1 (63) 



while the generating functional ()57p becomes (see also [54 



Q-f[A(p)|a(p),5^)] 



dp\{p) 



(64) 
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Taking functional derivatives of the logarithm of ()64p . the first, second and third order cumulants of 
the reference density are 

Ac-f (pi I S^) = p'-^f (pi I S^) = p{pi I S^) (65) 
k''°^(Pi,P2 I 5jv) = p'''Hpi,P2 I 5jv) - p{pi I Sn) p{p2 I Sf,) (66) 

- 1 ) p{pi I cS^) p(p2 I S^) (67) 

k""^ (Pl , P2 , P3 I 5jv ) = p'''^ (Pl , P2 , P3 I 5jv ) (68) 

;P2 I <Sm)p{P3 I '5]v) + 2/)(pi I Sm)p{P2 I <Sm)p{P3 I '5jv) 

3 - S^-f + 2 U(pi I 5^.)p(p2 I 5^.)p(p3 I 5^) (69) 

where the square bracket [3] indicates the number of distinct permutations which must be taken 
into account. The terms in the rounded brackets are readily recognised as the normalised factorial 
cumulants of the n distribution for a given fixed 





and so generalisation to arbitrary orders is immediate 



\"=0 / A=0 



l^'^'^^Pl, • • • ,Pr I 5jv) = KrN Yl P(.Pk I ^n)- (71) 

k=l 

This can be proven generally by defining the functional A[A(p)] = / dp\{p) p{p\Sm)/ {n)^ which has 
only a first nonzero functional derivative 5K/ 6\{pi) = p{p\Sf^) / {n)^ and the multiplicity generating 
function Z(A | 5^) = En^niv(l - ^T, in terms of which Q''^^[A] = ^[A[A]]. 

3.3 Internal cumulants for fixed 

Eq. (jTip shows that the differential cumulants of the reference distribution are directly proportional 
to the integrated cumulants K^-^n of n, which are zero only if TZun is poissonian. For fixed N , neither 
the integrated cumulants -fCr-,jv nor the differential ones are zero. While this has long been recognised 
in the literature |14| , the inevitable consequence was not drawn, namely that "poissonian" cumulants 
for fixed 

k{Pi,P2 I Sn) = p{pi,P2 I Sn) - p{pi I Sm) p{p2 I Sn), (72) 

k{Pi,P2,P3 I 5^r) = p{pi,P2,P3 I 5^r) 

- [3] p{pi ) P2 I '~>N ) P(P3 I <Sn) + 2p(pi I Sn) P{P2 I Sm) p{p3 I Sn) (73) 

etc. cannot possibly represent true correlations because they are nonzero even when the momenta are 
fully independent. It is known that the theory of cumulants needs improvement on a fundamental 
level which reaches well beyond the scope of this paper [68], [69], but those difficulties are irrelevant 
here. A first step which does address the above concerns was taken in Ref. [17J, where it was shown 
very generally on the basis of generating functionals that correlations for samples of fixed are best 
measured using the internal cumulants , which are defined as the differences between the measured 
and the reference cumulants of the same order 

K^{pi, ...,Pr\SN) = k{pi, ...,Pr\SN)- K'^'^iPl, ■■■,Pr\ Sn)- (74) 
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For our averaged-multinomial reference case, the internal cumulants of second and third order are 
given by the differences between Eqs. (j72]) and (j67|) and between ([73]) and (f69|) . resulting in 



H^{Pl,P2 I Sn) = p{pi,P2 I 5jv) - F2nP{Pi I Sm) P{P2 I 5jv), (75) 

+ Gsn p{Pi I Sn)p{p2 I Sn)p{p3 I ^at), (76) 



with 



LrSiV — •J-r2]v — -C^Sjv — -J —2 -—3 

and so on for higher orders. These internal cumulants are identically zero if and when the measured 
densities for fixed Sun are multinomials since then from Eq. ()62p 



p{pi,P2 I 5jv) p''''f(pi,p2 I ^at) = F2nP{Pi I ^at) p(P2 | ^at) (78) 

SO that k\pi,P2 I 5jv) — )• whenever the data is multinomial, while 

P{Pl,P2,P3 I Sn) Fsj^p{pi I Sn)p{P2 I 5Ar)/)(p3 I ^at) (79) 

ensures that {pi,p2,P3 \ S^) — )• in the same case. On another level, the internal cumulants always 
integrate to zero over the full good-azimuth space A, irrespective of the presence of correlations, 

/ dpidp2 k\pi,P2\Sn) = / dpidp2dp3 k^{pi,P2,P3\Sn) = (80) 
Ja ja 

and so on for all orders. Both properties will remain valid after transformation from three-momenta to 
invariant four-momentum differences in Section [3.41 In the case of Poissonian statistics, F^-n = 1 Vr, 
so that the above internal cumulants revert to their usual definitions. 

As stated in Section 13.11 the measured correlations may in addition be made insensitive to the 
one-particle distribution through normalisation. As set out in Ref. [T^, such normalisation is achieved 
for fixed A'^ by dividing the internal cumulants by the corresponding reference distribution density, 
which for the case at hand is given by Eq. (|62p . This leads to the second-order normalised internal 
cumulant 



K {pi,p2\Sf,) = r-^ - 1, 

F2N p{pi\Sm) p[P2 5]v) 



(81) 



while in third order we get 



r^I(„ ^ ^ \ q \ 1 f piPl,P2,P3\SM)-[3]p{Pl,P2\SM) piP3\SM) \ , ^F2M . . 

^ [pi,P2,P3 On) = 7 — rT-—? — rT~n — tfi + — ^- ^^-^^ 

F3N V P{Pi Sm)p{P2 Sm)p{P3 Sm) J FsiV 



3.4 Correlation integrals for momentum differences 

In femtoscopy, correlations are mostly expressed in terms of pair variables K = ^ {pi + P2) and differ- 
ence q = Pi —P2 or the invariant four-momentum [70] Q = \/—{pi — ^2)^ = \/ {Pi — P2Y — {Ei — £'2)^ 
where the energies are on-shell, E^ = \/Pr ~^ ■ As shown in Ref. ^52j, the formulation of eventwise 
counters as sums and products of Dirac delta functions makes it easy to change variables. Writing 
Prjv(pi, • • • ,Pr) as shorthand for p{pi, . . . ,Pr\ S^) etc., the second-order unnormalised internal cumu- 
lant in terms of Q is, for example, found from the identity / dQ kI^{Q) = f dQ J_^dpidp2 '«2Jv(PiiP2) 



15 



S{Q - ViPi-P2y-iEi-E2y) to be 

'^LiQ) = {yAQ-QtD) -F2.{{yAQ-Q'i')) ) m 





Nb'Na 

= P2n{Q) - F2N PlN<^PlNiQ), (84) 

where the counters in the second hne are defined by the terms in the first, while Q^j* = [—{P^ — Pj')'^]^^'^ 
and Q^j = [— (-Pf — Pj)'^]^^'^ are four-momentum differences between sibhng pairs aa and event mixing 
pairs ab respectively. It is easy to show that dQ P2n{Q) = ("'~)jv ^'^'^ dQ pii^®pii^{Q) = (n)^ 
and hence, as before, dQ kI^{Q) = 0, which will be true for any correlation whatsoever. The 
double event averages in the product term 

P.N^PiNiQ) = ((yS{Q-Qf^)\ \ (85) 

\ \ ''^ ' Nb'Na 

are the theoretical foundations of event mixing |52]; the inner 6-average is usually shortened to a 
smaller "moving average tail" subsample of ■ 

In third order, the "GHP average" invariant is defined as the average of three two-momentum 
differences over all pairs (with or without the V^), 

Qa = V-(.Pi - P2? - {P2 - P3? - {P3 - Pi? IV^ ; (86) 

it is related to the the invariant mass of three pions M3 = [pi +P2 +^3)^ by Q? = — m^. Other 
"topologies" such as the "GHP max" = Y^max[— (pi — P2?, —{P2 — P3?i —{ps — Pi?] can also be 
employed. For large multiplicities, the "Star" topology may be preferred [71j, but we shall not pursue 
it here. For the GHP average, the third internal cumulant is given by 

^L(QJ = {yHQ^- Q7^) ) - 3(( E E ^ - ^S) ) ) (87) 

'•^ ' Na Nb Na 

Nc Nb Na 



with Qf^^ = ^Jl[-{Pt - Pj? - {Pj - Pk? - {Pk - Pt?\^ and similarly for Ql^'^ and Q^^^j'. Second 
and third-order cumulants are normalised by, respectively, 

F2N P.N®pUQ) = F2MllY.^{Q-Qf-)) \ , (88) 

FSN PlN^PlN^PlNiQa) = F3N ///^ ^ " ) ) ) ' ^^^^ 

^'^^^ Nc Nb Na 

After transforming from momenta to Q, the formulae of Section [3.31 become 

«^2iv(<3) = P2n{Q) - F2nPiN^Pin{Q), (90) 
«^3Jv(Q«) = PsNiQa) - [3]p2N'S>PlNiQa) + G3N PlN® PlN® Pin{Q a) , (91) 

while the normalised cumulants of Section 13.31 become 



Ki{Q\SN)= ^ - 1, (92) 

F2nPin®Pin{Q) 



-Ifr^ I C \ - PiN{Qa) - [^]P2N®PlN{Qa) . „ -?^2jV 



3V^a| N) F^N PlN®PlN®PlN{Qa) F^n ^ ' 
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3.5 Effect of fixed-A^ correction factors 



To get a feeling for the size of the corrections involved, we measured the correction factors Fr^ and G^^f 
with the same UAl dataset and the same cuts as in Fig. 1. As shown in Fig. 2, the consequence of the 
clearly subpoissonian multiplicity distributions shown in Fig. 1 is that these factors are significantly 
less than 1, in contrast to the usual factorial moments of the charged multiplicity distribution which 
are superpoissonian with factors exceeding 1. For very low multiplicities N < 10, normalised internal 
cumulants are hence larger than their poissonian counterparts but converge to them with increasing A'^. 
Nevertheless up to ~ 30 corrections of more than 5% for K2 and more than 20% for compared 
to their uncorrected counterparts can be expected. By contrast, the additive correction G^n does 
not deviate much from the poissonian limit of 2 except for very small N. By contrast, unnormalised 
internal cumulants ()90p - ()9ip are far less sensitive to the multinomial correction. 




Figure 2: Left panel: Correction factors = {n{n—l))^/{n)^ (red squares) and F^^^ = {n(n—l){n—2))^/{n)^ 
(blue circles) as defined in Eq. (155)) as well as (green triangles) of Eq. ([TT]) . for the UAl dataset used in 
Rcf. [BT. Right panel: inverse factors as used in the normalisation of internal cumulants Eqs. (|8ip - ([5^ and 

dna-dMi). 

It is of interest to zoom in on the approach to the poissonian limit of 1 and to compare these 
corrections to the equivalent charged-multiplicity-based ones, which for the case of fixed A^, would be 
just N{N—1)/N'^ and N(N—l){N—2)/N^. In Fig. 3, the poissonian limit corresponds to zero on the 
y-axis. It is clear that the fixed- A^ factors go some way to correct for the fixed- A^ conditioning; the 
gap between them is approximately determined by {n-)^/N-, i.e. by the exact definition and outcome 
space for n. 



3.6 Eliminating fluctuations in n 

We return briefly to the first choice in Section [3.2.31 i.e. a{p \ Snn) = p{p \ SnN)/'^ which would permit 
different physics for each (A^, n) combination. If we were willing and able to do analyses for each SnN^ 
we would use the fixed-(A^, n) equivalent of (f75]) - ([76]) derived from K^'^^{pi^ . . . ,pr I Sun) = (— 1)''^^ {f — 
1)! ™nLiP(Pfc |5„^^), 

K^{Pl,P2 I Sun) = p{Pl,P2 \ Sun) " (l " ^) P{Pl I Sun) p{P2 \ Sun) (94) 
l^^{Pl,P2,P3 I Sun) = p{Pl,P2,P3 I Shn) - [3] p{pi,P2 \ Sun) p{P3 I Sun) 

+ 2{l- p{pi I Sun) P{P2 I Sun) p{P3 I 5„iv) (95) 

and normalise by p^'^^{pi, . . . ,pr \ Sun) = {n-/rf) nfc=i p{Pk \ Sun)- Where that is not possible, we 
can still average over the above to form "Averaged Internal" (AI) correlations (see Section [5]), but in 
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Figure 3: Comparison, for tlie same UAl data and cuts, of (1/-F2n) — 1 (red filled squares), (l/F^fj) — 1 (blue 
filled circles), both averages over n for fixed N, with prefactors [N'^ /N{N—1)] — 1 (red open squares) and 
[N^ /N{N—l){N—2)] — 1 (blue open circles). The poissonian limit is represented by zero on the y-axis, i.e. 100 
times the y-scale represents the percentage deviation from the poissonian limit. 



this case averaging over n for fixed A^, 

k'''°(Pi,P2 I S^) = {k^{pi,P2 I SnN))^ (96) 
K'^'"^(Pl,P2,f>3 I 5jv) = {k^{Pi,P2,P3 I ■Sun))^ (97) 



/JV 



and normalise if necessary by p^^'^{pi, . . . ,Pr\SM) = {{n-/ri'') YYk=iPiPk\SnN))f^- Given that this 
involves products of moments in the subsubsample 5„jv, event mixing would have to be restricted to 
the same subsubsamples also; for example 

N 

{{l-r)p{p,\Sn.)p{p2\Sn.)l = Y.^{l-l)j^ p{pi\n,N,P^)p{p2\n,N,P'') 

n=2 a,b£S„N 

(98) 

The transformation to pair variables works in the same way as in previous sections. 



4 Statistical errors 

While the various versions of internal cumulants, constructed above may all be relevant at some 
point, we concentrate on finding expressions for experimental standard errors for the unnormalised 
and normalised internal cumulants of Eqs. (j90P "(|93 p . This turns out to be more subtle than merely 
applying a generic root-mean-square prescription. We shall show in this section that standard errors 
implemented thus far may have been underestimated even in the standard two-particle case. 

The calculations performed in this section belong to the "frequentist" view of probability; a proper 
Bayesian analysis, which can be expected to rest on more solid foundations, is beyond the scope of 
this paper. The two viewpoints can reasonably be expected to yield similar results in the limit of large 
bin contents and sample sizes. 

In this section, we often simplify notation by writing pr{Q \ S^) — )• PrN and pi0pi{Q \ S^) — )• pfj^ 
etc, since the formulae apply to samples and variables of any kind. 
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4.1 Propagation of statistical errors 



Because cumulants can be measured only through the moments that enter their definitions, the first 
task is to identify which moment variances and covariances are needed. By means of standard error 
propagation ^68j, we find the sample variances for second-order internal cumulants of Eqs. (j9U|) and 
Eq. ([92]) to be 



vav{Ki{Q I 5jv)) = var(p2jv - -^2jv pIn) = var(p2jv) + var(/>^jy) - 2 cov(p2jv, /oD (99) 
var(i^2^(Q I Sn)) = var - l] (100) 

\F2N PtN ) 



1 P2N 

F^N Pin 



var(/j2iv) , var(pf 



P2N 



+ 



PfN 



, COv(p2iV,/)D 
2 

P2N ■ PiJV 



(101) 



under the assumption that var(i<2jv) is much smaller than the other variances, so that F2N can be 
treated as a constant; this is the case if there are many bins for Q. Similarly, from ()76p . the variance 
of the unnormalised internal cumulant is, assuming G3 of Eq. (|77p to be constant. 



var(K^((5, I Sn)) = var [p^N - ^P2nPin + G3 pf^^] 

= var(/>3iv) + 9var(/92ivPijv) + G|var(/3f^) 
— 6 cov (p.jN, P2NP1N ) + 2G3 COv(p3jv, /?ijv 



6G3 COV {p2NPlN,Pir 



while the normalised version has variance (again assuming var(G3) <^ vaipr 



var(K|((3„ I Sn)) = var 

PSN 



PiN — 3/32JvPlJV + G^p^f 



FsNPtf 



/?2 



var 



3p2ivPi 



FsnP^ 



PSN 



IN 

3p2ivPl JV 



var(p3jv - 3p2NPi 
{p3N - 2>p 

2nPiN J 



+ 



PsN 3/?2JvPi 

var(/9^^) 

FlJV 



+ G3 



(102) 
(103) 

(104) 



2cov(p3jv - 3p 

2N PiNi Pin) 
{psN - 3/? Pin)p IN 



Fsn Pfr. 



+ 



var(p3jv) + 9var(p 

2nPiN 

) - 6cov{p3n,P2nPin) _^ vav{pf^) 

{psN - 3/? 
6cov {p2N Pin, pIn) - 2cov(p3jv,/)^^) 

{PsN ~ 3P2nPin)Pin 



(105) 



The unnormalised cumulants ([90|) - ([9T]) and their variances (j99|) and ()103p require knowledge of the 
multiplicity factorial moments F2N , -Psjv , so that the individual terms must be accumulated until the 
entire sample has been analysed. By contrast, the normalised cumulants (j92p - (j93p and their variances 
(jlOip and ()105p contain the multiplicity moments only as prefactors. 



4.2 Expectation values of counters 

While we shall not make direct use of the results in this section, it is nevertheless useful briefly to 
consider what we might mean by an "expectation value of experimental counters and densities" . For 
any scalar function f{p) of the momenta, the theoretical expectation value E[f] is defined as the 
integral over the entire outcome space of / weighted by a "parent distribution" P{p), an abstract 
entity supposedly containing everything there is to know on this level, 

E[f{p)] = I dpP{p)f{p). (106) 
Jn 

Purely theoretical concepts such as P{p) and E[f] should be given little or no room in a strongly 
experimentally-oriented study. In calculating standard errors on counters below, we shall, however. 
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make use of the exact factorisation that expectation values provide whenever two variables x, y are 
statistically independent, E[xy] = f dx dy P{x,y) xy = E[x]E[y]. 

Expectation values for pairwise variables such as the four-momentum difference Q we are consid- 
ering here must be based on the underlying physics. We can deduce some properties of the parent 
distribution based on the usual definition of the femtoscopic correlation function 

is a function of two entirely different quantities: the four- momentum differences Q'*" of "sibling" 
tracks taken from the same event a and one constructed from the mixed-event sample using tracks 
from different events, written as Q*^^, Q^'^ etc. For second-order correlations, the parent distribution is 
therefore necessarily a two-variable probabilitjlll P(Q"'*,Q^^) which, depending on whether the cases 
h=a and c=a occur, may or may not factorise into a "sibling" and a "mixed" marginal probability 

P(Q-,Q^^) = P,(Q-)P^(Q^-) iffa/6/c, (108) 
but (unless a=h=c) the marginals will always be 

P,{Q^^) = J dQ^'P{Q'"',Q^''), (109) 

Pra{Q^l = j dQ'"'P{Q'"',Q^''). (110) 

The shapes of Ps{Q) and Pm{Q) must necessarily be different since it is precisely this difference that 
leads to a nontrivial signal in (|107p . In terms of this joint probability, we can write expectation values 
of eventwise counters (separately for inclusive, fixed- or fixed-n cases) as 



(111) 



^W^)]= / dQ-dQ'^P(Q-,Q''^)p,, = ^ / dQ'^-PM''')m'^--Q%) = Y,PM%). (112) 
Later, we shall meet expectation values for cases such as a=c, 



E[p{Q''^) /)(Q'^^)] = Y.J2 ^^"^ ^^"^ J^(Q"", Q"') 5(Q'^'^ - gf;)5(Q'^' - Q 



ab\ 
kl) 



which definitely does not factorise. The above expressions can be simplified because we know that the 
parent distribution is not a function of the individual track indices i,j,k,i 

P{Qt^,QM) = P{Q^'',Q'l yhj,k,e (114) 

and similarly Ps(Q?/) = i^s(<5"") and PmiQ^) = Pm{Q''^)- For the event -averaged counters, this 
results in 

E[paam = {n^lPsiQn (115) 

EiPbciQ)] = {n,l{n,lPM'') (116) 



We could argue that there are three different variables Q"", Q" and Q where the last two differ in the sense that 
Q"** contains a track from the "current" event while Q*"^ does not. As shown below, this distinction is unnecessary as 
long as we keep careful track of possible occurrences of equal event indices. 
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or in terms of the notation of Section 13.41 

E[p,^m = {n\Ps{Qn (117) 
E[p m = {n)^Pm{Q"'). (118) 

As mentioned, we do not need the factorisation (jlOSp of P{Q°'°', Q^^) as long as we keep careful track 
of the equal-event-indices cases. Whenever a^h or ot^c, independence of the events ensures that 
expectation values of products of any functions f{Q"'"') and g{Q°'^) of the pair variables do factorise, 

E[f{Qn giQ'^')] = ElfiQ'^'^)] i^b(Q'^')] a / b (119) 

For third-order correlations, the parent distribution is a function of three different variables Q""", 
Qbbc g^j^j^ Q'^^f containing respectively three, two or one track from the same event and corresponding 
considerations regarding equal and unequal event indices apply there, too. 



4.3 Statistical error calculation from first principles 



It was shown in Section 14.21 that expectation values would have well-defined meanings in terms of 
underlying parent distributions and their marginals if their parent distributions were known, which, 
however, they are not. We are therefore forced to revert from expectation values E[-\ to sample 
averages (•) after completing a calculation. The real use of such expectation values in frequentist 
statistics has been in the form of a gedankenexperiment which we now reproduce from Kendall |68j . 
Let X be any generic eventwise counter or any other eventwise statistic. Since the formulae in this 
section remain true for inclusive and fixed-A^ samples, we omit any notation related to N in this 
derivation. In this simplified notation, the well-known standard error of the sample mean is given 
by (simplifying £ — 1 ^ £) 



a 



(120) 



which follows from the combinatorics of equal and unequal event indices by the above artificial use of 
expectation values, reverting from expectation values E[-] to sample means (•) in the last step: 



var((x» = E[{xy] - E[{x)]^ = ^2 E " ^t^"^] 

a,b 

= ^2^2 [^l^aXb] - E[Xa] E[xb]\ +^2^2 [^[^J ^[^b] - ^[^'^l ^[^b] 
a=b a^b 



(121) 



1 r 



(122) 



where we have used the fact that E[xaXi,] = E[xa] E[xi,] Wa^b and assumed that all x are identically 
distributed, E[xa] = E[x] Va. Equality or inequality of event indices is thus crucial. We shall follow the 
same approach below, keeping careful track of equal and unequal event indices, factorising expectation 
values for unequal event indices, and reverting to sample means in the last step. 



4.4 Variances and covariances for multiple event averages 
4.4.1 Statistical errors for second-order cumulants 

According to Eqs. (f99|l - (|105p . we must handle variances and covariances of products of several event 
averages. To derive these, we shall use the following shortened notation: Letting 6f^ = 6{Q — Qfj) etc. 
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then paa = Tl^i^j ^{Q~Qi^) = ^ij eventwise pair counter for event a, while pbc = X^i j ^iQ~ 

Qij) — Ylij^ij is the mixed-event counter of events b and c (with by^Cy^a assumed), so that P2{Q) = 
{Paa}^ while pi0pi{Q) = {pab)^i^ is a double event average. We reserve the event index a for the 
"sibling" event whose correlations are currently being analysed, and use indices b,c, . . . ,t,u,v,w, . . . 
for events entering the event-mixing partsU All quantities are assumed to be measured within a 
particular subsample but we omit the A^-subscript and the argument. The event-index subscripts 
such as ( )bc above are included or omitted depending on whether they convey relevant information on 
the specific averaging. 

In this notation, the method that led to Eq. ()122p reads 



var(p2) = var{{paa)J 



a.b 



Pbb) - E{paa) E{p 



bb) 



[{Pa 



{PaaY 



(123) 



The same method of disentangling the combinatorics of equal and unequal event indices is applied 
consistently to all variances and covariances below. The pf-tevm in second-order cumulants has vari- 
ance 

1 



var{pi0pi) 



var 



E{pbc)E{pde) 



(124) 



b^c d^e 

The case b^^c^d^e yields zero, but the cases b=d^c^e and three other equivalent combinations yield 



var(pi(g)pi, 



(£:2)2 L 
A£'^ 



E[pbcPbe] - 
{pbcPbe)bce 



E[pbc] E[pb, 

- {pbc) {Pbe) 



{PbcPbe) - {Pl®Pl) 



(125) 



) and in the thirqj assumed £ ^ 1. Note that the 
{{T.i^j^\^)c{T.k^i^M)e)b simplified to the 



where in the second step we reverted E[\ - 
requirement b^c^e implies that {pbcPbe)b^^ 

square of a single counter {{Yli=/=j ^ij)'l)b' event mixing involves three different events, not two. 
Secondly, the combinations of two equalities b=d^c=e and b=e^c=d in (jl24p yield another term of 
order £~'^, 



2 

^2 L 



{pbdpbd)bd - ipl^plf 



(126) 



which we can safely neglect when {v?'^/£ <C 1 except when there are few bins or large multiplicities 
even in small bins. It is worth emphasising that the extra factor 4 which appears in Eq. (|125p arises 
from the same method that has been used for decades to justify use of Eq. (|120p . We find, by the 
same method, that the covariance between p2 and pi®pi is given by 

cav{p2.,pi®pi) = cov{{pdd), (pbc)) = 7772 [E{pddPbc) - E{pdd) E{pbc) 



££1^^\ 

d b^tc 



2 r 
£ 



PddPdc 



{p2){pl(^Pl) 



{Pdd){Pdc)\ 

SO that we must in addition accumulate, for every event d, the product of the counters 

i^j k 



(127) 
(128) 



dc\ 

kl I ■ 

I c 



(129) 



^While it is irrelevant whether event a is included or excluded in theoretical calculations of event mixing, it should 
never be used in actual implementations of mixing. 

^Due to the factorisation of the expectation values earlier on, the fact that index b appears in two separate sample 
averages does not prevent us from replacing {ptc) and {pbe) by (/9i®pi)'^- 
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Note that there is no restriction on track indices k^i or k^j in the d-event, meaning that events with 
n(a) = 2 contribute to this counter which would otherwise not be the case. Combining these, we find, 
to leading order m. and renaming mixed-event indices, 



Var(4) = i |((/5aa - '2.F2NPac) {Paa " Pad)) ^ " (P2 " '2,F2n P2{Pl® Pljf'^ 



(130) 



with all event indices strictly unequal and c- and d-event averages understood where appropriate^ 
Contrasting this with the traditional way to calculate the same variance, 

Var(K^) = ^ |((/Oaa - F2NPac) {Paa " F2MPad)) " {P2 " F2N P2{Pl® Pljf^ (131) 

it is clear that in previous analyses the two possible ways to set a equal to 6 or c were overlooked, 
while normal (non-internal) cumulants also omit the -F2jv- 

4.4.2 Statistical errors for third-order cumulants 

In third order, we shall need 6'^^^ = 6{Qa — Qijj!.) and similar quantities, and the notation for counters 
Paaa, Paab and Pabc Corresponding to the event averages psiQa) = (paaa), P2®Pi{Qa) = {Paab) and 
Pi<8'/0i<8)/9i(Qa) = {pabc) respectively. Clearly, ay^by^c must hold in the third order case. We obtain for 
the necessary third-order quantities (shuffling and/or renaming indices if necessary) 



Var(p3) = ^[{{Paaa)\ - pI 
(.P3,P2^Pl) = ^ ^ \^{prrTPsst) " ^(/Srrr) F{hst) 



r sj^t 



^[{prrrPrrs) + (^PrrrPrss 



(132) 



(133) 



with PrrrPrrs = ^i^j^k^Vik^i^mi^n^lm-n) ^ ^'^^ PrrrPrss = Hi^j^k^ljk^f.{^m^n^lmn') ^- ^he 

remaining variances and covariances needed for third-order correlations with GHP topology are, after 
renaming of indices. 



r 

Var(/92«>Pl) = ^ ^ \E{pggePccd) - E{pgge) E{pccd) 



gyte c^d 

■A r 

^ [{pggdip 



while we neglect 



1 



'ggc + Pgcc + Pddc + Pdc 



PggcPggc + PggcPgcc - '^{P2®Pl] 



4(p2®Pl)^ 



The next term is simpler, 

'(P3, P?) = ^ ^ \E{pmpuvw) - E{pttt) E{p^ 



GOV 



t UJ^V^W 



3 

£1 



[PtttPtuv) - P3Pl 



(134) 



(135) 



(136) 



While this factorised form is instructive, it cannot be used directly since F2N can be determined only on completion 
of the entire sample analysis. Each of the counter products in (|130p must hence be implemented separately. 
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but the following is not, 



cov 



{p2'^Pl,pf) = -^f^ ^ ^ ^E{puuvPxyz) - E{puuv) E{pxyz 
Uy^v xj^yj^z 

3 ^{PwwxPwyz) + {PwxxPwyz) — '^{P2Pl){p\) 



3S^ 



+ 



(137) 

Pwxy) - {p2pl){pl)\, 

(138) 



and the large number of combinations makes the variance of p\ particularly complicated, 



var(pi(X'/9i(g)pi) 



+ 



_6_ 



:3^ X] X] \_^iPdehPbpq) - E{pdeh) E{pbpq 

18£:^ - 



(139) 



■\{pbehPbpq)- {Pl){p\) 



+ 



{pbehPbeq) - {pl){Pl 



{pbehPbeh)- {p\){pl) 



(140) 



For large 8, the leading order terms will usually dominate, so that we can neglect the subleading 
termsO To leading order, we therefore obtain after substitution in ()102p and again omitting brackets 
for non-a event averages 



Var(K|) = ^1 ^Paaa + ^Paab{Paac + Pace + Pbbc + Pbcc) + 'dC'lpaaaiPabcPade) 

— QPaaaiPaab + Pabb) + QG^PaaaPabc — '^^G^PaabiPacd + Pbcd 

I a 

pI + 36{p20pif + 9Gl (p?)2 - 12p3ip2^pi) + 6G3 psipif - 36Gl{p2^pi)ipif 



(141) 



{Paaa — '^Paab — "^Pabb + '3G^Pabc){Paaa — '^Paad — "^Padd + "^G^Pade, 



P3 - Qp2®Pl + 3G^{pi] 



(142) 



While the factorised form is again instructive, it cannot be calculated in this form within the a-loop 
in the analysis since the G3JV constants are known only on completion of the entire sample analysis. 
Rather, the full palette of product counters pp has to be accumulated and averaged and combined 
only in the final phase of the analysis. 

5 Averaged internal cumulants 

As is only an approximation for the true total event multiplicity anyway, and for cases of small 
sample statistics, it may be necessary or desirable to group subsamples of fixed into multiplicity 
classes A^ S -B]. It is important, however, not to simply lump all events within this multiplicity 
class into a single "half-inclusive" subsample, because, as has long been known [14] , that results in 
terms entering the cumulants which arise solely to "multiplicity mixing" (MM) of events of different 
A^. Given the arbitrary choice of such MM correlations are spurious and avoided in favour 



If and when large bins are used and the sixth power of the measured positive-pion muhiphcity becomes comparable 
to £, subleading terms will have to be included. This requirement is less trivial than it may sound, since for subsamples 
of fixed multiplicity A'', the number of events En is much smaller than £, while of course n may be substantial when A'^ 
is large. For UAl, £m = 0(10") while £ = O(10'*). 
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of "Averaged-Internal" (AI) correlation^ defined as follows. Using the renormalised multiplicity 
distribution 

7^; = = (143) 

the AI unnormalised cumulants, reference distributions and normalised AI cumulants are 

B 

4\Q \Sab)=Y. ^'iv ^2(<3 I Sm), (144) 

N=A 
B 

Pi®Pi{Q I 5ab) = K F2N Pi®Pi{Q I 5jv), (145) 

i^2^^(Q I Sab) = "^'^^''^f\ = J:^n=aT^'^P2{Q\S.) ^ 

Pi®Pi{Q\Sab) En=aT^'mF2nPi^Pi{Q\S^) 

Note that the correction factors Fon, which are normalised factorial moments of n for fixed A^, are 
part of the summed normalisations!^ In third order, we have correspondingly 

B 

\Sab)=Y1 ^'n 4 (a \Sn), (147) 
N=A 
B 

Pi^Pi^PiiQa I Sab) = ^ 'R-'n Fsn Pi(^Pi®Pi{Qa I S^), (148) 

N=A 

Ki^iQ. I Sab) = f^^^j^^l Y (^^9) 

Pl®Pl®Pl{Qa\SAB) 

Note also that ()148p holds for the normalisation only and not for the last term in Kg, which is 
(3F2iv— -Fsiv) Pi®Pi®Pi — but that is already taken care of in the formula ([76]) for Kg itself. Expressions 
for an inclusive (all-A) multiplicity summation of internal cumulants are obtained from the above by 
setting ^4=0 and B = oo. The AI (Averaged Internal) correlations Eqs. ()144p and ()147p represent 
refined versions of what has traditionally been termed "Short-Range Correlations" , differing from the 
original formulae [14] [16] by the i^2jv and Gsjv factors respectively. This was originally pointed out in 
Ref. [17] but only for multinomials in A^. 

Regarding variances and standard errors for AI correlations, we first note that, since subsamples 
5jv are strictly mutually independent, a variance over the [A.B] range is simply the weighted sum of 



^Historically, this issue was discussed under the name "Short- Range Correlations" and "Long-Range Correlations" 
[14) . Since current usage of the term "Short-Range Correlations" refers to correlations over small scales in momentum 
space, we rather define them more accurately as "Averaged Internal" (AI) correlations and "Multiplicity-Mixing" (MM) 
correlations, noting also that the correction factors in Eqs. (|144|) - (|149|l do not appear in the earlier literature. 
An expression with a correction factor outside the sums such as 

EiV^«^2(Q|5„) 



(n(n-l)) Ejv^'ivPi®Pi(<3!'5iv) 

is inconsistent with AI correlation averaging if the single-particle spectra or some other physical effect change significantly 
within the range [^4,5]. We also note that the above differs from the formula used in Ref. t64j for second-order correlations 
in q. In the present notation, the cumulant used in Ref. |64) reads 

Eiv^'ivP2(q|5N) 



Eiv 7^'„^Pl®/^l(q|<S„) 

i.e. a correction for an A'^-multinomial rather than the weighted sum of n-multinomials used in Eqs. (|144p - (ll46p . 
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the corresponding fixed- variances. From Eqs. ()144p and ()149p we have to ah orders r 



B 



i^fiQ I Sab) = Yl ^'^ (<3 I ^iv) r = 2, 3, . . . , (150) 



N=A 



and given the independence of any functions / and g of different multiplicity subsamples, E[f{SN) 
= E[f{Si^)] ■ E[g{S^>)] V TV 7^ TV', we conclude that 



var[Kf{Q I Sab)] = ^iK)^ var[4iQ I 5^)] (151) 

AT 

varK(Q I 5ab)] = J](7^;)2 (F,^)^ var[p^(Q | 5^)] (152) 

N 

cov[k;:^^(Q I 5ab),/9UQ I Sab)] = Y^iK)^ cov[k^(Q | S^),p\{Q \ S^)] (153) 

N 

which are known functions in terms of Sections 14.4.11 and 14.4.21 while for the normalised cumulants in 
[A, -B], covariances between numerator and denominator are (omitting the Q), 

cov[«;2(5jv),Pi(5iv)] = coy[p2{Sm),Pi®Pi{Sn)] - F2ivvar(/)i(g)pi(5iv)) (154) 
co\[kI{Si,),p\{Sm)] = coy[p^{Sm),p\{Sm)] - 3cov[p2(»Pi(5jv),/0?(5jv)] + G3var[pf (5;v)] (155) 

so that the normalised range cumulants have variances 

y^t[K^\Q\Sab)] = {K^')'' 
and standard errors are given b;'^^' 



Yax[Kr{SAB)] ^ ysx[pi{Sab)] _ ^ GOV [k^ {Sab), Pi {Sab ) ] 
(uriSAB))'^ {Pi {Sab ) ) ^ i^r (Sab ) ■ Pi {Sab ) 



(156) 



a{K^' {Sab)) = ^w&t{K^'{Q \ Sab))- (157) 
6 Event mixing algorithms 

"Event mixing" [72j is widely used to simulate uncorrelated or semi-correlated quantities such as p\ 
and P2®pi- The idea has always been to use the experimental sample at hand to simulate the baseline 
of independence referred to in Section [3. II in such a way that criteria 2 (independence of momenta), 4 
(reproducing the one-particle momentum space distribution) and 6 (normalisation) are all addressed 
simultaneously. Ideally, all effects bar the desired correlation are elegantly removed in this way. 

For the internal cumulants and their variances and covariances derived above, event mixing requires 
keeping track of counters of all orders in each of the subsamples Sj^. A count of event indices in 
Section H] shows that in a brute-force calculation we would need, for each subsample, a minimum of 



^^One might expect g^K^^Sab)) to include a prefactor of the sort seen in Eq. (|120|l i.e. something like 
^J\va.x{K^' [Q I Sab)]/[B — A], but this would be incorrect. The reason in that the formulae H151|) - ()153[l for range 
AB can be considered as an average, so that we can apply the methods of Section [4.31 to obtain the same results. For 
example, considering H2' {Q\Sab) = /ti of (|144p as an average and writing kKQISn) = K2n, the variance on this 
average is 

var[Ki] = E[{-Ktf] - ^[Ki]^ 

= Eiv=iv'(^'«)' {E[fiL]-E[4{Q\S^)f) = Eiv(7^'«)'var(K2«) 
which is identical with (|151[) . Therefore, division by B — A is incorrect. 
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five independent event averages or 0{£^) event combinations; furthermore, caution would advise not 
to use the same event in calculating related counters, so that selection and use of more than five 
events in mixing is advisable. The resulting excessive number of full event averages, mixing every 
(sub)sample event with every other one, is therefore not feasible. 

If the order of events in the sample is random, the multiple event averages can be simplified by 
the use of the following multiple event buffer algorithm. 

1. A single overall event loop equivalent to the event index a runs over the entire inclusive sample 
S. A given event a will have a multiplicity N=N{a), so for that particular a correlation analysis 
for subsample Sj^ is advanced by one event while the others remain dormant El In this way, a, 
which always refers to the sibling event, runs over all events of every subsample S^. There 
is no need to either explicitly sort the inclusive sample into subsamples or to run multiple jobs 
for fixed A^. 

2. The first Sg event^^ of a given multiplicity are used solely to fill up the buffer without 
doing any analysis. Once a given buffer has been filled, event mixing analysis proceeds for that 
subsample as follows: 

(a) An newly-read a-event is assigned to the N=N{a) buffer, the earliest event in that buffer 
is discarded, and sums for averages entering and F^j^ as well as the sibling counters 
Paa, Paa, {Paaf and {paaaf are Updated. 

(b) Event combinations for mixed-event counters are built up by picking any one of the £"3 — 1 
other events in that buffer and calling it 6, thereafter picking any one of the remaining 
£"5 — 2 events in the buffer, calling it c and so on. While to third order only five events 
(including the sibling event) are needed to construct all the counters required, in practice 
it is better to use different mixing events for different counters to root out even traces of 
unwanted correlations between different mixing counters. The random selection of events 
rather than tracks for mixing is necessary to ensure that more than one track per event can 
be used as required for counters of Sections [3H1] such as phbc- 

(c) For a given event set b,c,d,..., the mixing counters are incremented using all possible 
combinations of the tracks in event a together with all the nb,nc,... tracks in the 
selected events 6, c, . . . mixing events. For example pbbc would use all possible nb{nh — 1) 
pairs of 6-tracka^^l together with all possible Uc single c-tracks. The mixing of all tracks of a 
given event rather than just selected ones ensures that the fluctuations in n for given fixed 
A^ are automatically contained in the counters. 

(d) For constant a, the process of selecting events 6, c, . . . is repeated Cmix = 10-100 times to 
reduce the statistical errors, avoiding events that have been used in previous selections. 
Efficiency can be maximised by tuning of both the number of events £b stored in each 
buffer and by the number of resamples Cmix- 

3. Once the entire sample has been processed via the a-event loop, the b- c- d-event averages are 
normalised by Cmixi^N — Sb), while the primary a- aver age is normalised by {S^ — £b)- 

4. Unnormalised and normalised correlation quantities, their standard errors and correction factors 
are constructed by appropriate combinations of averaged counters. 

'^^Inevitably, there are very few events in the high-multiphcity tail of the entire sample. These must be treated 
separately, e.g. by putting all events with N greater than some threshold into a single buffer. 
'^^The number of events in a buffer £b is usually kept the same for each buffer. 
^^As in the definition of the counters, each pair is counted twice: these are ordered pairs. 
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5. Results from fixed- A'^ subsamples can then be combined into AT Correlations over partial ranges 
of A*" or the entire inclusive sample at the end of the event loop using the methods outlined in 
Section [5i 



Discussion and Conclusions 

1. Correlations are only defined properly if the null case or reference distribution is defined on the 
same level of sophistication as the correlation itself. Translated into statistics, the six criteria 
set out in Section [3. II for a reference sample for correlations at fixed charged multiplicity N lead 
straight to the definition of the reference sample as the average of multinomials given in Eq. ()57p , 
weighted by the conditional multiplicity distribution TZnN- Assigning Bernoulli probabilities 
a{p I Sun) = {p{p I ^niv))^^/ {n)j^ yields the reference density ([62]) and generating functional 

2. From the theorem that internal cumulants are given by the difference between measured and 
reference cumulants we obtain normalised and unnormalised internal cumulants which satisfy 
every stated criterion for proper correlations for fixed-iV samples. 

3. We have highlighted the distinction between n, the particles entering the correlation analysis it- 
self, and N, the particles determining the event selection criterion for a particular semi-inclusive 
subsample. Various correction factors are shown to be fair to good approximations of these 
exact results in some cases but far off the mark in others. Surprisingly, normalised cumulants 
are far more sensitive to these corrections, through the normalisation prefactor, than their un- 
normalised counterparts. To belabour the point: For any variables (xi,X2), different definitions 
for correction factor 1 /F for fixed- correlations of positive pions. 



can be very important at low multiplicities, with F = 1 (Poisson) being the worst approximation, 
F = N{N—1)/N'^ a fair one, and {n{n—l))^/{n)^ being the best. 



Refs. \74\ [75] in an approach based on probabilities rather than densities. Ref. |76j specifically 
calls the inclusion of these correction factors meaningless because the theory then requires that 
the emission function be identically zero. We note that the argument in all those references relates 
to inclusive samples, while for the samples of fixed considered in this paper the prefactor, 
which is an average of n at fixed A^, is a necessity. Either way, the arbitrariness of the use or 
non-use of the prefactor has been eliminated here based solely on considerations related to the 
reference distribution. 

4. The problem posed in this paper, viz. the relation between charged multiplicity on the one hand 
and correlations based on the conditional n-distribution T^ruv, has attracted little attention in 
the literature. Indeed, almost all theoretical work on multipion correlations, as for example 
summarised in Ref. [77], starts from the projection of final-state events, with all their different 
particle species, onto the single-species subspace of either n+ (positive pions) or n_ (negative 
pions) correlations, to the exclusion of the other charge. It would be interesting to see a combined 
theory for both rij^ and n_, which would encompass all the work done so far plus correlations 
between unlike-sign pions and, of course, the issue raised by us here. 

5. As shown in Figs. 2-3, correction factors for third order are larger than the second-order ones. 
For higher r-th order correlations, the effect of using a fixed- A^ subsample is suppressed by 




(158) 



For inclusive correlations, correction factors such as (ra)^/(n(n— l))inci were proposed early on in 
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approximately l/(n^~^)jv for unnormalised cumulants but actually worsens for normalised cu- 
mulants due to the normalisation prefactors (n)^/(n-)jv for small n. The importance of accurate 
correction of normalised quantities therefore rises with order of correlation. 

6. The difference between poissonian and internal cumulants is largest at small multiplicities n. The 
mixed-multinomial prescription will therefore be required for any correlation analysis involving 
small n, independently of the magnitude of N. Apart from the usual suspects of leptonic, 
hadronic and low-energy collisions, the low-n case occurs both for very restricted phase space 
(such as in spectrometer experiments) and for correlations of rare particles such as kaons and 
baryons, even for large N. 

7. Since n fluctuates according to the conditional multiplicity distribution TZnN , the degree to which 
fixed- correlations differ from inclusive ones is strongly coupled to the character of TZnN- In 
general, TZnN is subpoissonian and so F^jv falls below the poisson limit of 1. The correction from 
fixed-A poissonian to internal normalised cumulants is hence upward, not downward as in the 
case of multiplicity-mixing corrections. 

8. While not the main subject of the present paper, some light is cast on the relationship be- 
tween three levels of correlation, viz. correlations inherent in the overall multiplicity distribution, 
multiplicity-mixing correlations, and the true internal correlations for fixed A^. Each of these 
can and should be treated separately. The Averaged- Internal correlations of Section [5] are a 
compromise solution which may be useful both for physics reasons and for small datasets. 

9. The fixed-A corrections discussed here are separate and complementary to other important 
effects at low multiplicity. Refs. [771 EH] highlight, for example, possible effects of "residual 
correlations" resulting from projecting from multipion to two-pion correlations. 

Energy-momentum conservation would also play a role. Borghini [79\ [80] has, for example, 
calculated the effect of momentum conservation for normalised two- and three-particle cumulants 
in momenta and pt- However, the saddlepoint method used applies to the large- A^ limit and 
the results cannot be directly applied to the low-A (and hence low-n) samples under discussion 
here. Indeed, momentum conservation will be near-irrelevant for cases of large A and small 
n as discussed above, but the small-n corrections of this paper will remain important. For 
the specific case of like-sign pion femtoscopy, the fact that only n ~ AN/2Q out of the A" 
charged pions are used and that momentum conservation constraints include all other final-state 
particles both imply that momentum conservation constraints may be less important than the 
internal-cumulant correction introduced here. 

For the specific choice of correlation variables Q and Qa for two- and three-particle cumulants, 
the contribution of momentum conservation to cumulants at small Q will be small since the 
counts will be dominated by pairs at small (At/), Ay) and intermediate {pti,Pt2)- As pointed out 
by [T], momentum conservation exerts greatest influence at large pair or triplet momenta and 
hence mostly at large Q, where it may lead to moments and cumulants which do not converge to 
a constant as presupposed in most fits. The ad hoc method of multiplying fit parametrisations 
by a prefactor 1 + cQ with c a free parameter does not adequately address the problem. 

Regarding the multiplicity dependence of the influence of energy-momentum conservation, Ref. [2] 
calculates the effect of energy-momentum conservation on single-particle differential observables, 
and finds significant systematic effects. No doubt this must also be the case for multiparticle 
observables, although, as we have pointed out above, the effect of conservation laws will be di- 
luted by the fact that fewer than half of the final-state particles of any given event are actually 
used in the present analysis. Detailed investigations are beyond the scope of this paper. 
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10. We have also recalculated statistical errors for products of event averages starting from the 
original prescription which forms the basis of frequentist statistical error calculations. Compared 
to conventional statistical error calculations, new prefactors appear in our calculations — see 
e.g. Eq. (jl25p and the results in Section [4.4.21 — which have surprisingly been missed so far. 

11. We note that the present formalism is still in the frequentist statistics mindset, which may be 
inaccurate for low multiplicities and should be supplanted by a proper Bayesian analysis. The 
final word has certainly not been spoken about correlation analysis of small-n datasets. 
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